/*
 *  fftw.cpp
 *  phaseit3d
 *
 *  Created by Anton Barty on 26/07/11.
 *  Copyright 2011 Anton Barty. All rights reserved.
 *
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include "phaseit3d.h"

cFFTWdata::cFFTWdata(){
	nx = 0;
	ny = 0;
	nz = 0;
	nn = 0;
	data = NULL;
}


cFFTWdata::cFFTWdata(long n){
	nx = n;
	ny = n;
	nz = n;
	nn = n*n*n;
	free(data); data = NULL;

	data = (fftwf_complex *) fftwf_malloc(nn*sizeof(fftwf_complex));
	if(!data) {
		printf("ERROR in fftwData: Malloc failed\n");
		exit(1); 
	}
	
}


cFFTWdata::~cFFTWdata(){
	
	fftwf_free(data); data = NULL;
	fftwf_destroy_plan(plan);
	fftwf_destroy_plan(iplan);
	fftwf_cleanup_threads();

}	



void cFFTWdata::create(long n){
	
	create(n,n,n);
	
}

void cFFTWdata::create(long nnx, long nny, long nnz){
	nx = nnx;
	ny = nny;
	nz = nnz;
	nn = nx*ny*nz;
	free(data); data = NULL;

	data = (fftwf_complex *) fftwf_malloc(nn*sizeof(fftwf_complex));
	if(!data) {
		printf("ERROR in fftwData: Malloc failed\n");
		exit(1); 
	}
}



/*
 *	Convert 3D (xyz) triple into 1D array index
 *	x is fastest axis, z is slowest axis
 */
long cFFTWdata::element(long x, long y, long z){	
	return x + nx*y + (nx*ny)*z;
}



/*
 *	Create single threaded FFTW plan
 */
void cFFTWdata::createFFTWplan(void) {
	
	createFFTWplan(1);

}

	
/*
 *	Create multi-threaded FFTW plan
 */
void cFFTWdata::createFFTWplan(int	nthreads) {

	char	wisdomFile[1024];	
	FILE	*fp;

	/*	
	 *	Sanity check
	 */
	if(nx<=0 || ny<=0 || nz<=0) {
		printf("Illogical FFT dimensions: %li x %li x %li\n", nx,ny,nz);
		exit(1);
	}
	
	
	/* 
	 * 	Threads
	 */
	printf("\tSetting nthreads = %i\n",nthreads);
	if (fftwf_init_threads() == 0) {
		printf("\t\tCould not initialise threads\n");
		exit(1);
	}
	fftwf_plan_with_nthreads(nthreads);
	
	
	/* 
	 *	Read Wisdom from file 
	 *	Size of fftwf_complex used to determine whether we are using fftwf_ or fftwf_ 
	 */
	strcpy(wisdomFile,getenv("HOME"));
	if (sizeof(FFTW_DATA_TYPE) == 2*sizeof(float)) 
		strcat(wisdomFile,"/fftwf_wisdom");
	else if (sizeof(FFTW_DATA_TYPE) == 2*sizeof(double)) 
		strcat(wisdomFile,"/fftw_wisdom");
	
	printf("\tImporting FFTW wisdom from %s\n",wisdomFile); 
	fp = fopen(wisdomFile, "r");
	if (fp != NULL) {
		if (!fftwf_import_wisdom_from_file(fp) )
			printf("\t\tError reading wisdom!\n");
		fclose(fp); 	/* be sure to close the file! */
	}
	else {
		printf("\t\tCould not open FFTW wisdom file %s\n",wisdomFile);
	}
	
	
	/* 
	 *	Generate FFTW plans 
	 */ 
	printf("\tCreating forwards plan....\n");
	plan = fftwf_plan_dft_3d(nx, ny, nz, data, data, 1, FFTW_PATIENT);
	printf("\tCreating inverse plan....\n");
	iplan = fftwf_plan_dft_3d(nx, ny, nz, data, data, -1, FFTW_PATIENT);
	
	
	
	/* 
	 *	Export wisdom to file 
	 */
	fp = fopen(wisdomFile, "w");
	if (fp != NULL) {
		printf("\tExporting accumulated wisdom to file %s\n",wisdomFile);
		fftwf_export_wisdom_to_file(fp);
	}
	fclose(fp);
	
	
}



/*
 *	Do the FFT
 */
void cFFTWdata::fft(int direction){
	if(direction == 1)
		fftwf_execute(plan); 
	else if (direction == -1)
		fftwf_execute(iplan); 
	else {
		printf("Error: Undefined FFTW direction\n");
		exit(1);
	}
	
}


/*
 *	Speed test
 */
void cFFTWdata::speedTest(int nit){
	clock_t			start, start2;		
    time_t			startt, endt;
    float			dt;
    float		avg_t;

	
	printf("\nCalculating timing for %i FFT/iFFT pairs\n",nit);
	start = clock();
	time(&startt);
	
	for(long it=1; it<=nit; it++) {
		index_loop {
			data[p][0] = 1.0;
			data[p][1] = 0.0;
		}
		
		start2 = clock();
		fftwf_execute(plan); 
		fftwf_execute(iplan); 
		
		time(&endt);
		dt = endt-startt;			
		//dt = difftime(endt,startt);
		printf("\t%li : %li^3 FFT/iFFT pair took %3.2f sec\n",it,nx, dt/it);
		
		avg_t += (float)(clock()-start2);
		//printf("\t%i : %i^3 FFT/iFFT pair took %3.2f sec\n",it,nn, (clock()-start2)/(float)CLOCKS_PER_SEC);
	}
	
	time(&endt);
	dt = endt-startt;			
	//dt = difftime(endt,startt);
	avg_t = (clock()-start2)/ nit;
	avg_t = avg_t / ((float)CLOCKS_PER_SEC);
	printf("%li^3 FFT/iFFT pair average (CPU time): %3.3f sec\n",nx, avg_t);
	printf("%li^3 FFT/iFFT pair average (human time): %3.3f sec\n",nx, dt/(float) nit);
	
	
}


